rm(list=ls())

setwd("C:/Users/bschneer/Documents/GOV1/GOV 2001/Replication/attributes")
library(cem)
library(mvtnorm)
library(foreign)
library(apsrtable)
library(Matching)
sims<-500000
set.seed<-02138

cf<-read.dta("prepare/camp_finance_clean.dta", convert.factors=F)
names(cf)<-sub("_",".",names(cf))

#Campaign Finance CEM

match<-cf[,c("judge.vote","treatment","l.treatment","r.treatment","circuit","dec.year","jcs","year.birth","minority.judge","gender.judge","lower.dir","abs.jcs","jud.experience")]
match["circuit"]<-as.factor(match[,"circuit"])

imb1<-imbalance(match$treatment,match,drop=c("judge.vote","treatment","l.treatment","r.treatment"))

dec.year_cut<-c(1980,1991.5)
jcs_cut<-c(-0.63, -.3, 0,0.3)
year.birth_cut<-c(1910,1929.5)
minority.judge_cut<-c(0,1)
gender.judge_cut<-c(0,.5)
lower.dir_cut<-c(0,.5)
abs.jcs_cut<-c(0, .5)
jud.experience_cut<-c(0,1)
circuit_gps<-list("1","2","3","4","5","6","7","8","9","10","11","DC")

cp<-list(dec.year=dec.year_cut,jcs=jcs_cut,year.birth=year.birth_cut,gender.judge=gender.judge_cut,minority.judge=minority.judge_cut,lower.dir=lower.dir_cut,abs.jcs=abs.jcs_cut,jud.experience=jud.experience_cut)
gps<-list(circuit=circuit_gps)

cem.match<-cem(treatment="treatment",data=match,cutpoints=cp,grouping=gps,drop=c("judge.vote","l.treatment","r.treatment"))
cem.match
imb2<-imbalance(match$treatment,match,drop=c("judge.vote","treatment","l.treatment","r.treatment"),weights=cem.match$w)

m1<-MatchBalance(treatment~jcs + year.birth + gender.judge +  minority.judge + lower.dir + abs.jcs + dec.year, data=match)
m2<-MatchBalance(treatment~jcs + year.birth + gender.judge +  minority.judge + lower.dir + abs.jcs + dec.year, data=match, weights=cem.match$w)

before1<-c()
after1<-c()
for (i in 1:7) before1<-cbind(before1,m1$BeforeMatching[[i]][3:6])
for (i in 1:7) after1<-cbind(after1,m2$BeforeMatching[[i]][3:6])

before1<-data.frame(before1)
after1<-data.frame(after1)

names(before1)<-c("jcs","year.birth","gender.judge","minority.judge","lower.dir","abs.jcs","dec.year")
names(after1)<-c("jcs","year.birth","gender.judge","minority.judge","lower.dir","abs.jcs","dec.year")
before1;after1
save(before1,after1,file="cfmeans.RData")

#diagnostic<-imbspace(cem.match,match,M=250,depth=3,plot=F)

#pdf("cf_match.pdf",width=6,height=6)
#plot(diagnostic,data=match,explore=F)
#dev.off()

spec1<-formula(judge.vote ~ treatment)
spec2<-formula(judge.vote ~ treatment + jcs + year.birth + gender.judge + minority.judge + lower.dir + abs.jcs + as.factor(dec.year))
#Drop gender.judge in spec2a due to inflated variance
spec2a<-formula(judge.vote ~ treatment + jcs + year.birth + lower.dir)
spec3<-formula(judge.vote ~ l.treatment + r.treatment)
spec4<-formula(judge.vote ~ l.treatment + r.treatment + jcs + year.birth + gender.judge + minority.judge +lower.dir + abs.jcs + as.factor(dec.year))

model1<-glm(spec1,binomial(logit), data=match)
model2<-glm(spec2,binomial(logit), data=match)
model3<-glm(spec1,binomial(logit), data=match,weights=cem.match$w)
model4<-glm(spec2a,binomial(logit), data=match,weights=cem.match$w)

model5<-glm(spec3,binomial(logit),data=match)
model6<-glm(spec4,binomial(logit), data=match)

#Hypothesis Test
v.l<-summary(model6)$cov.unscaled[2,2]; v.r<-summary(model6)$cov.unscaled[3,3]; cv.lr<-summary(model6)$cov.unscaled[2,3]
beta.l<-summary(model6)$coefficients[2,1]; beta.r<-summary(model6)$coefficients[2,2]
print(abs(beta.l-beta.r)/(sqrt(v.l+v.r+2*cv.lr)))

apsrtable(model1,model2,model3,model4,stars=0,model.names=c("Full: Naive","Full: Multivariate",
"Matched: Naive", "Matched: Multivariate"),notes=list("se.note"),omitcoef=9:27,coef.names=c("(Intercept)","Treatment","Judge Ideology","Year of Birth","Female Judge","Minority Judge","Lower Court Direction","Circuit Ideology"))
apsrtable(model5,model6,stars=0,model.names=c("Full: Naive","Full: Multivariate"),
notes=list("se.note"),omitcoef=10:28,coef.names=c("(Intercept)","Treatment (L)","Treatment (R)","Judge Ideology","Year of Birth",
"Female Judge","Minority Judge","Lower Court Direction","Circuit Ideology"))

#"Some covariates (Gender and Race) omitted from Matched model due to low frequency of outcomes"

invlogit <- function(x){
  1 / (1 + exp(-x))
}

ate1 <- function(model, sims){
  data        <- model$data
  beta        <- coef(model)
  V           <- summary(model)$cov.scaled
  beta.draw   <- rmvnorm(sims, mean=beta, sigma=V)
  pr.men      <- matrix(0, nrow(as.matrix(data)), sims)
  pr.women    <- matrix(0, nrow(data), sims)
  X.men       <- X.women <- model.matrix(model)
  X.men[,2]   <- 0
  X.women[,2] <- 1

  diff.holder <- matrix(0, nrow(data), sims)

  for(i in 1:nrow(data)){
    p0   <- invlogit(X.men[i,] %*% t(beta.draw))
    p1   <- invlogit(X.women[i,] %*% t(beta.draw))
    diff <- p1 - p0
    
    diff.holder[i,] <- diff
  }
  return(apply(diff.holder, 2, mean))
}

e<-matrix(0,nrow=4,ncol=3)
for (i in 1:4){
ate<-ate1(get(paste("model",i,sep="")),1000)
e[i,]<-cbind(mean(ate),quantile(ate,.025),quantile(ate,.975))
}

data1<-data.frame("Campaign Finance",c("Full Naive","Full Multivariate","Matched Naive", "Matched Multivariate"),e)
names(data1)<-c("Issue Area","Type","ATE","Lower CI","Upper CI")

#Capital Punishment CEM

cp<-read.dta("prepare/cap_pun_clean.dta", convert.factors=F)
names(cp)<-sub("_",".",names(cp))

match<-cp[,c("judge.vote","treatment","l.treatment","r.treatment","circuit","dec.year","jcs","year.birth","minority.judge","gender.judge","lower.dir","abs.jcs","jud.experience")]
match["circuit"]<-as.factor(match[,"circuit"])

imb1<-imbalance(match$treatment,match,drop=c("judge.vote","treatment","l.treatment","r.treatment"))

dec.year_cut<-c(1980,1990,2000)
jcs_cut<-c(-0.63, -0.22, 0.18, 0.58)
year.birth_cut<-c(1910, 1918.8, 1927.6, 1936.4, 1945.2, 1954)
minority.judge_cut<-c(0,.5)
gender.judge_cut<-c(0,1)
lower.dir_cut<-c(0,.5)
abs.jcs_cut<-c(0, 0.87)
jud.experience_cut<-c(0,1)
circuit_gps<-list("1","2","3","4","5","6","7","8","9","10","11","DC")

cp<-list(dec.year=dec.year_cut,jcs=jcs_cut,year.birth=year.birth_cut,minority.judge=minority.judge_cut,gender.judge=gender.judge_cut,lower.dir=lower.dir_cut,abs.jcs=abs.jcs_cut,jud.experience=jud.experience_cut)
gps<-list(circuit=circuit_gps)

cem.match<-cem(treatment="treatment",data=match,cutpoints=cp,grouping=gps,drop=c("judge.vote","l.treatment","r.treatment"))
cem.match
imb2<-imbalance(match$treatment,match,drop=c("judge.vote","treatment","l.treatment","r.treatment"),weights=cem.match$w)

m1<-MatchBalance(treatment~jcs + year.birth + gender.judge +  minority.judge + lower.dir + abs.jcs + dec.year, data=match)
m2<-MatchBalance(treatment~jcs + year.birth + gender.judge +  minority.judge + lower.dir + abs.jcs + dec.year, data=match, weights=cem.match$w)

before2<-c()
after2<-c()
for (i in 1:7) before2<-cbind(before2,m1$BeforeMatching[[i]][3:6])
for (i in 1:7) after2<-cbind(after2,m2$BeforeMatching[[i]][3:6])

before2<-data.frame(before2)
after2<-data.frame(after2)

names(before2)<-c("jcs","year.birth","gender.judge","minority.judge","lower.dir","abs.jcs","dec.year")
names(after2)<-c("jcs","year.birth","gender.judge","minority.judge","lower.dir","abs.jcs","dec.year")
save(before2,after2,file="cpmeans.RData")

#diagnostic<-imbspace(cem.match,match,M=300,depth=3,plot=F)

#pdf("cp_match.pdf",width=6,height=6)
#plot(diagnostic,data=match,explore=F)
#dev.off()

spec1<-formula(judge.vote ~ treatment)
spec2<-formula(judge.vote ~ treatment + jcs + year.birth + minority.judge + gender.judge + lower.dir + abs.jcs + as.factor(dec.year))
spec3<-formula(judge.vote ~ l.treatment + r.treatment)
spec4<-formula(judge.vote ~ l.treatment + r.treatment + jcs + year.birth + minority.judge + gender.judge + lower.dir + abs.jcs + as.factor(dec.year) )

model1<-glm(spec1,binomial(logit), data=match)
model2<-glm(spec2,binomial(logit), data=match)
model3<-glm(spec1,binomial(logit), data=match,weights=cem.match$w)
model4<-glm(spec2,binomial(logit), data=match,weights=cem.match$w)

model5<-glm(spec3,binomial(logit),data=match)
model6<-glm(spec4,binomial(logit), data=match)

apsrtable(model1,model2,model3,model4,stars=0,omitcoef=8:26,model.names=c("Full: Naive","Full: Multivariate","Matched: Naive", "Matched: Multivariate"),notes=list("se.note"),
coef.names=c("(Intercept)","Treatment","Judge Ideology","Year of Birth","Minority Judge","Female Judge","Lower Court Direction"))
apsrtable(model5,model6,stars=0,model.names=c("Full: Naive","Full: Multivariate"),notes=list("se.note"),omitcoef=10:16,coef.names=c("(Intercept)","Treatment (L)","Treatment (R)",
"Judge Ideology","Year of Birth","Female Judge","Minority Judge","Lower Court Direction","Circuit Ideology"))

#Hypothesis Test
v.l<-summary(model6)$cov.unscaled[2,2]; v.r<-summary(model6)$cov.unscaled[3,3]; cv.lr<-summary(model6)$cov.unscaled[2,3]
beta.l<-summary(model6)$coefficients[2,1]; beta.r<-summary(model6)$coefficients[2,2]
print(abs(beta.l-beta.r)/(sqrt(v.l+v.r+2*cv.lr)))

e<-matrix(0,nrow=4,ncol=3)
for (i in 1:4){
ate<-ate1(get(paste("model",i,sep="")),1000)
e[i,]<-cbind(mean(ate),quantile(ate,.025),quantile(ate,.975))
}

data2<-data.frame("Capital Punishment",c("Full Naive","Full Multivariate","Matched Naive", "Matched Multivariate"),e)
names(data2)<-c("Issue Area","Type","ATE","Lower CI","Upper CI")

rule <- c(4,3,2,1)

dplotl <- function(data){
outside <- (data[,4] > 0 | data[,5]<0)
	plot(data[,3], rule, ylim=c(0.5,4.25),xlim=c(-.4,.4),xlab="", ylab="",yaxt="n", col="black", lwd=3, pch=20 )
abline(h=1, lty=3, col="lightgray")
abline(h=2, lty=3, col="lightgray")
abline(h=3, lty=3, col="lightgray")
abline(h=4, lty=3, col="lightgray")
abline(v=0, lty=1)
segments(data[,4],rule,data[,5],rule, col="black", cex=1.5, lwd=2)
segments(data[,4][outside],rule[outside],data[,5][outside],rule[outside], col="red", cex=1.5, lwd=2)
points(data[,3][outside], rule[outside], col="red", lwd=3, pch=20 )
axis(side=2, at=rule, labels=data[,2], las=1, cex=0.5, lwd=0)
	}

dplotm <- function(data){
outside <- (data[,4] > 0 | data[,5]<0)
	plot(data[,3], rule, ylim=c(0.5,4.25),xlim=c(-.4,.4),xlab="", ylab="",yaxt="n", col="black", lwd=3, pch=20 )
abline(h=1, lty=3, col="lightgray")
abline(h=2, lty=3, col="lightgray")
abline(h=3, lty=3, col="lightgray")
abline(h=4, lty=3, col="lightgray")
abline(v=0, lty=1)
segments(data[,4],rule,data[,5],rule, col="black", cex=1.5, lwd=2)
segments(data[,4][outside],rule[outside],data[,5][outside],rule[outside], col="red", cex=1.5, lwd=2)
points(data[,3][outside], rule[outside], col="red", lwd=3, pch=20 )
	} 

pdf(file="judge_characteristics.pdf", width=7, height=4)
par(mfrow=c(1,2),oma=c(0,5,0,0))
#CF
dplotl(data1)
title("Campaign Finance")
mtext("Full N: 117; Matched N: 34",cex=.75)

#CP
dplotm(data2)
title("Capital Punishment")
mtext("Full N: 375; Matched N: 155",cex=.75)
dev.off()

#Also Create Indidvidual Plots

#CF
pdf(file="cf_judge_characteristics.pdf", width=7.25, height=6.5)
par(mfrow=c(1,1),oma=c(0,5,0,0))
dplotl(data1)
title("Campaign Finance")
mtext("Full N: 117; Matched N: 34",cex=.75)
dev.off()

#CP
pdf(file="cp_judge_characteristics.pdf", width=7.25, height=6.5)
par(mfrow=c(1,1),oma=c(0,5,0,0))
dplotl(data2)
title("Capital Punishment")
mtext("Full N: 375; Matched N: 155",cex=.75)
dev.off()

